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Abstract 

We estimate the probabihty of detecting a gravitational wave signal from coalescing compact 
binaries in simulated data from a ground-based interferometer detector of gravitational radiation 
using Bayesian model selection. The simulated waveform of the chirp signal is assumed to be a 
spin-less Post-Newtonian (PN) waveform of a given expansion order, while the searching template 
is assumed to be either of the same Post-Newtonian family as the simulated signal or one level 
below its Post-Newtonian expansion order. Within the Bayesian framework, and by applying a 
reversible jump Markov chain Monte Carlo simulation algorithm, we compare PN1.5 vs. PN2.0 and 
PN3.0 vs. PN3.5 wave forms by deriving the detection probabilities, the statistical uncertainties 
due to noise as a function of the SNR, and the posterior distributions of the parameters. Our 
analysis indicates that the detection probabilities are not compromised when simplified models are 
used for the comparison, while the accuracies in the determination of the parameters characterizing 
these signals can be significantly worsened, no matter what the considered Post-Newtonian order 
expansion comparison is. 
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I. INTRODUCTION 



Kilometer-size ground-based interferometric detectors o 



ravitational radiation have be- 
From locations in 



i Q, 3, Q. 



come operational at several laboratories around the world 
the United States of America, Italy, Germany, and Japan, these instruments have started 
to search, in the kilohertz frequency band, for gravitational waves emitted by astrophysical 
sources such as spinning neutron stars, supernovae, and coalescing binary systems. 

Among the various sources of gravitational radiation that these instruments will attempt 
to observe, coalescing binary systems containing neutron stars and/or black holes are ex- 
pected to be the first to be detected and studied. These signals have a unique signature 

n 

that enables them to be extracted from wide-band data by digital filtering techniques [5|. 
This signature is their accelerating sweep upwards in frequency as the binary orbit decays 
because of energy loss due to the emission of gravitational radiation. Coalescing binaries 
have a potential advantage over other sources in signal-to-noise ratio (SNR) by a factor that 
depends on the square-root of the ratio between the corresponding number of cycles in the 
wave trains ^. 

The standard technique used for extracting these "chirps" from the noisy data is called 
Matched Filtering. In the presence of colored noise, represented by a random process n{t), 
the noise-weighted inner product between the data stream d(t) recorded by the detector and 
the gravitational waveform template h{t) is defined as 

2 Re (1, 



where the symbol "over d and h denotes their Fourier transform, S{f) is the one-sided power 
spectral density of the noise n(t), fi, fu are the limits of the frequency band of interest, 
and the * represents complex conjugation. From this definition the expression of the SNR 
can be written in the following form Q 

By analyzing the statistics of the SNR^, it is possible to make statements about the presence 
(or absence) of such a gravitational wave signal in the data. This operation of course needs to 
be performed over the entire bank of templates over which the SNR statistics is built upon, 
since a gravitational wave signal is in principle determined by a (finite) set of continuous 
parameters. 



The effectiveness of tlie matclied filtering procedure relies on the assumption of exactly 
knowing the analytic form of the signal (possibly) present in the data. Recent breakthroughs 
in numerical relativity 

Baa 

have started to provide a complete description of the radiation 
emitted during the inspiral, merger, and ring-down phases of generic black hole merger 
scenarios. Although the ability of obtaining numerically all the templates needed in a data 
analysis search (perhaps hundred of thousands of them) might be practically impossible, 
in principle we should be able to compare these numerically derived waveforms against 
various analytic templates obtained under different approximating assumptions. Work in 



this direction has already started to appear in the literature [id, llll, within the so called 
"frequentist framework" , in which estimates in the reduction in SNRs and inaccuracies of the 
determination of the parameters characterizing the signal, due to the use of approximated 
waveforms, have been derived. Depending on the magnitude of these degradations one can 
decide whether to use these approximated waveforms as templates in a data analysis search. 

Since it can be argued that contiguous PN approximations should well characterize the 
differences between the "true signal" present in the data and the highest-order PN approx- 
imation, in this paper we perform such a comparison within the Bayesian framework. An 
analogous, frequentist analysis has recently been performed by Cutler and Vallisneri [13] for 
the case of super massive black holes binaries observed by LISA (the Laser Interferometer 
Space Antenna). Their approach relied on the use of the Fisher-Information matrix, which 
is known to give good results in the case of large (hundreds to thousands) SNRs. In the case 
of ground-based interferometers instead, since the expected SNRs will be probably smaller 
than 10, a parameter estimation error analysis based on the Fisher-Information matrix would 
lead to erroneous results jl^ . 

In this paper we will estimate the loss in probability of detection (i.e. loss of evidence 
of a signal to be present in the data) as a function of the SNR in the following two cases: 
(i) the true signal present in the data is a spin-less PN3.5 waveform and the search model 
is represented by a spin-less PN3.0, and (ii) the true signal is a spin-less PN2.0 waveform 
and the model is a spin-less PN1.5 waveform. We have limited our analysis to these two 
separate cases in order to cover the region of the PN approximations that have already, or 
are in the process of, being used in the analysis of the data collected by presently operated 
ground-based interferometers. 

Our approach relies on a Bayesian Markov Chain Monte Carlo (MCMC) technique, as 
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MCMC methods have successfully been applied to a large number of problems involving 
parameter estimations jl5| in experimental data sets. In our analysis the chirp signals (the 
one present in the data and the one used as the model) are characterized by five parameters: 
the two masses of the system, mi and m2, their time to coalescence tc, the coalescence phase 
0(7, and their distance r from Earth. 

Bayesian MCMC methods have already been proven to be capable of estimating the five 
parameters of a PN2.0 chirp signal embedded into noisy data of a single interferometer 
when using a PN2.0 based model 16|, and in a coherent search in the time domain for nine 
parameters using a PN2.5 model [17| and a PN3.5 model ISj for the phase. However, it 
has never been shown before how the resulting evidence and probability distributions of 
the parameters are affected by the usage of different PN models. Here we will estimate 
the evidence of a signal being buried in noise and, at the same time, derive the probability 
distribution of its parameters when the gravitational wave form of the model is a simplified 
version of the signal present in the data. A Bayesian analysis naturally justifies Occam's 
Razor Q due to the penalization of unreasonably complex models by the integration 
over the parameter space. The paper is organized as follows. In Section [III we provide 
a brief summary of the Bayesian framework and its implementation in our problem. After 
deriving the expressions of the likelihood function and the priors for the parameters searched 
for, in Sec. IIIII we describe the Markov-Chain Monte Carlo sampling technique adopted for 
calculating the posterior distributions. We then specify in Sec. |IV]the different simulations 
we have conducted by introducing the wave forms, noise specifications, and parameter sets. 
The final results of our simulations are presented in Sec. |Vl displaying the MCMC based 
posterior distributions for the parameters and involved models. The estimated posterior 
probabilities of the models are presented as a function of SNR with the corresponding 
uncertainties over the noise realizations. We find that the difference in detection probability 
when using a simplified model rather than the true one is negligibly small in comparison to 
these uncertainties. The posterior credibility regions of the parameters reveal offsets from the 
true parameter values that can be much larger than the statistical uncertainty, for both the 
PNl. 5/2.0 and PN3.0/3.5 model comparisons. The PN2. 0/2.0 and PN3.5/3.5 comparisons 
on the other hand, always yield credibility regions that cover the true parameter values. 
Finally, in Sec. IVII we provide our comments and concluding remarks. 
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II. THE BAYESIAN FORMULATION 



In this section we will derive the Bayesian full probability model for our problem, which 
involves the comparison between the two possibilities of having either a signal and noise 
or just noise in the data. This requires the determination of the likelihoods, the prior 
distributions for the parameters associated with the models, and the resulting posterior 
distributions. 

A. Model definition 

Let us suppose we observe a data stream d{t) — St{0; t) +n{t) containing the instrumental 
noise n{t) and a chirp signal St{d; t) that we will regard as the "true^^ signal. Here, 6 is the 
vector representing all the parameters associated with the signal, and the noise is assumed to 
be a stationary Gaussian random process of zero mean. In the Fourier domain the observed 
data can equivalently be written as d{0; f) = St{6; /) +n(/) (where tilde denotes the Fourier 
transform operation) and we will refer to this expression as model A4t- 

In what follows we will assume the "true" signal St{d; f) to be the gravitational wave 
emitted by a coalescing binary system and represented by a spin-less Post-Newtonian ap- 
proximation in phase and Newtonian in amplitude for which — {mi, m2, r, to, (pc}^- Here 
mi and m2 are the masses of the rotating objects, tc is the coalescence time, r the absolute 
distance to the binary system, and (f)c the phase of the signal at coalescence. 

We will then describe the detection and estimation of the parameters of the "true" signal 
by relying on a spin- less lower-order Post- Newtonian waveform, Ss{d; /). This simpler model 
will be referred to as model Ms- 

The derivation of the detection probability implies a comparison between model ^-nd 
the null-model, which postulates mere noise h{f) within the data. This model will bereferred 
to as model Ain- Note that no parameter enters into this model. 

B. The likelihood function 

Since we have assumed the distribution of the random process associated with the noise 
of the detector to be Gaussian of zero-mean, it follows that the likelihood function is pro- 
portional to the exponential of the integral of the squared and normalized residuals between 
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d{Ot] f) and signal template s{6; f) over the frequency band of interest [fiyfu]- We de- 
fine the data d{f) := d{6t; f) to be modeled by Ait with its ''true" parameter vector Ot. 
The comparison with the simplified model Aig, results in the following expression for the 
likelihood function 

pid\Ms, 6) oc exp (^-2 ^^^^^^lyf^d/ j , (3) 

where Sn{f) is the one-sided power spectral density of the noise. 

By substituting d{f) = St{f, Ot) + n{f) into Eq. [3l the likelihood function becomes 

mM., 0) oc exp (-2 £ mM±m^MSl,f) . (4) 



s,.(f) 

In analogy to Eq. |H under model Ain (with no parameters) the likelihood assumes the 
following form 

p{d\Mn) oc exp (^-2 -^-jj-^ df j . (5) 

For comparison reasons, in the case of using the ''true" model the likelihood function 
becomes 

'■^^ m,Ot)+n{f)-St{f,0)\' 



p{d\Mu9)^exp\^-2j^ ' '''' ' g'J^j^ '''' " d/J , (6) 

which will then give information about the impairment in detection when using model M.s 
instead of M.f Note, that the investigation of this question requires to do the two different 



model comparisons separately, i.e. vs. M.s and M.n vs. A^j [21| 



The next step needed for completing our Bayesian full probability model is the identifi- 
cation of suitable prior distributions for the five parameters characterizing the chirp signal. 
Our derivation will closely follow that described in {igI . \\\ . 



C. Prior distributions 

The derivation of appropriate priors p{0) bears significant influence on the evidence of a 
signal presence within noise, since the prior identifies the size of the parameter space which 
the evidence is based on. 

A detailed description of the derivation of the prior distributions, and in particular for 



the masses mi, m2, and distance r, can be found in 16|, ll7|]. In short, the masses are 



assumed to be uniformly distributed over a specified range, [mminii^max]) and the prior 



distribution for the distance is chosen to be a cumulative distribution of having systems out 
to a distance x smaller than r, P{x < r), proportional to the cube of the distance, x^. In 
order to obtain a proper prior distribution that does not diverge once integrated to infinity, 
it is down-weighted by including an exponential decaying. This accounts for the Malmquist 
effect 22| and includes the assumption of uniform distribution for the masses. The resulting 



distribution function p{mi,m2,r) can be written as follows 
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p(mi,m2,r) OC /[m_„,m™„,](mi) /[m_„,m^„,](m2) 



1 + exp I 1254^1 V (7) 



where 



^ / (mi+m2)V6r ' 
and I[mmin,nT'max]{m) is the so-called "Indicator function", equal to 1 when rrimin < m < 
TTLmax and zcro elsewhere. The latter term containing the log-amplitude in the sigmoid 
function of Eq. [7] is the down weighing term mentioned earlier, and it depends on two 
constants a and h. These are determined by requiring a smooth transition of a mi,m2 
insp iral system being detectable with two specified probabilities at two given distances. In 
[m \\\ a and h are determined by choosing a (2 — 2)Mq inspiral system to be detectable 
with probabilities 0.1 and 0.9 out to distances 95 Mpc and 90 Mpc respectively. In our work 
we will also make such a choice. 

Fig. [1] shows the joint prior distribution of the masses mi, m2, and marginal distribution 
of the distance r, using Eq. [3 Although initially a uniform distribution is assigned to the 
masses, the conjunction of distance and masses results in a higher detectability of large 
masses. The number of possible binary systems increases quadratically with the distance 
but the down weighing of the prior is significantly seen above 500 Mpc as we allow masses of 



up to 5OM0. In 16|, ll7|] for example, the masses where restricted to SMq and therefore the 



prior values on the distance were much smaller with a distribution mode at around 75 Mpc. 



As far as the time to coalescence is concerned, we have assumed it to be uniformly 
distributed over a time interval of 1 second centered around the value identified by the 
masses of the binary and the lower frequency cut-off of the detector [23|. This search 

n 

range for the time to coalescence tc is larger than that used in pj| because when using 
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FIG. 1: Joint prior distribution of the masses mi, m2, and marginal distribution of the distance r, 
using Eq. [71 Although initially a uniform distribution is assigned to the masses, the conjunction of 
distance and masses results in a higher detectability of large masses. The number of possible binary 
systems increases quadratically with the distance. The down weighing of the prior is significantly 
seen at around 500 Mpc and above as we allow masses of up to 50Mq. 

the simplified model the posterior peak can be offset from the true value by more than the 
posterior width. 

Finally, we chose the phase of the signal at coalescence, (pc, to be uniformly distributed 
over the interval [0,27r], i.e. p(0c) = I[o,27r[(0c)/(27r). 

The choice of priors is different when it comes to the analysis of model Ain- Since 
this model postulates mere noise, there are no parameters entering the likelihood, which is 
therefore a constant. 

The final remaining step in defining the Bayesian procedure is to assign prior probabilities 
to the models themselves. Since we have no a priori knowledge, the unbiased choice is equal 
probability for each. 
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D. Posterior distribution 



By applying Bayes' theorem using the hkehhoods and priors defined above, we then derive 
the multidimensional posterior probability distribution for the model and its parameters 



p{i, Oi\d) 



p{Mn) ■ v{d\Mn) if^ = 

p{Ms,e)-p{d\Ms,e) if 2 = 1 



(9) 



jp{Ms^ 6) ■ p{d\Ms^ O)d0 + p{Mn) ■ pid\Mn) 

where i G {0, 1} corresponds to the two models {Ain, -Ms}- In the same way, it is possible 
to derive the posterior for the comparison of Ain vs. Ait- A Bayesian analysis naturally 



justifies Occam's Razor 19|, l20l] due to the penalization of unreasonably complex models by 



integrating over the parameter space resulting in the preference for a simpler model. 



III. BAYESIAN MODEL SELECTION AND PARAMETER ESTIMATION 



There exist various techniques for tackling this multi-dimensional problem. One possible 



approach is to calculate the so called Bayes factors 



21 



24 



25| . which are the ratios of 



the glo 
(BIG) 



Da 



likelihoods of the models that are involved. The Bayesian Information Criterion 
26| can be used as an approximation to the Bayes factor. However, it is possible to 
address the problem of sampling from the multidimensional posterior distribution in Eq. [9] 
by implementing a relatively new procedure, called Reversible Jump MCMC (RJMCMC) 
technique 27|, l28|, which simultaneously addresses the problems of model selection and 
parameter estimation. The RJMCMC is combined with traditional fixed dimension MCMC 
techniques that sample from the parameters of the current model. In the following we will 
briefly review the MCMC algorithm that we will use in our analysis. 



A. Metropolis Coupled Markov chain Monte Carlo 



In Simulated Tempering 



2^ . the "temperature" becomes a dynamic variable on which 



a random walk is conducted during the entire sampling process. The joint distribution of 
temperature and remaining parameters, however, requires the normalization constants of 
the distributions given the temperature. Other approaches like the Tempered Transition 
method [s^ or the Metropolis-Coupled chain (a.k.a. parallel tempering algorithm) jsil do 
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not need normalization constants. The latter approach has been advocated in the astro- 



physical literature 



21| and it has also been implemented in 17|, |32 1 



In a Metropolis-Coupled chain [31i], sampling is done in parallel from k different distri- 
butions pj{i,6i\d),j G {1, . . . , A;}. The real posterior distribution of interest is denoted by 
Pj{i,6i\d) with parameter vector 6i, whereas the distributions of higher orders j > 1 are 
chosen in such way that the sampling process is facihtated. Usually, different temperature 
coefficients are applied [sj] that flatten out the posterior modes. During the sampling from 
the k distributions, from time to time, attempts are made to swap the states of a randomly 
chosen pair of distributions. 

The posterior in the present context can be regarded as a canonical distribution 

p{Mn) ■ exp (-2/3, Md/) if z = 



where 



+ p(X„) ■ p(X„) ■ exp (-213, J'" ^■§jJJ<iA (11) 



with inverse temperatures (3j,j G {!,..., k}. For higher values of j, the posterior modes are 
flattened out and the sampling process is eased. A temperature scheme for our Metropolis- 
Coupled chain uses A; = 10 different Pj values with j G jj^-^^,k}, where /3i = 1 is the 



temperature of the original posterior distribution. As in [17|, [3^ , the prior distribution is 
purposely not involved in the temperature scheme as the prior information at high tempera- 
tures is preserved. Eg. [TO] converges to the prior distribution if /? ^ whereas a temperature 
scheme, had it been applied to the entire posterior distribution, would merely yield a uniform 
distribution. 

The inverse temperatures l3j,j G {1, . . . ,k} are unknown parameters that must be deter- 
mined prior to each simulation. It is obvious that the highest temperature needs to account 
for the nature of the likelihood surface. The stronger the signal, the higher the modes and 
the more likely it becomes for the MCMC sampler to get trapped. The acceptance prob- 
ability of a proposed jump in a basic Metropolis-Hastings algorithm is determined by the 
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product of the ratios between the proposals, the priors, and the hkehhoods [15|, 13^, l35|] of 
the proposed parameter vector and current state parameter vector. For a coarse assessment 
of the nature of the posterior surface, we can neglect the prior distribution as it is much 
smoother than the likelihood surface. With symmetric proposals, the likelihood ratio is 
therefore the key factor in analyzing the depth of the modes in the posterior surface. 

Although in the simplified model comparison the parameter estimates can be far off the 
true parameter values, in the true model comparison the true parameter values are expected 
to be good estimates of the parameters. This fact can be used for a coarse assessment of the 
nature of the likelihood surface. Since the log-likelihood of the true model at the true pa- 
rameter values is logLHt = —2 f^^ ^^^^^^^^^^^^^df, it can be compared to the log-likelihood of 
the null-mo del logLHn = -2 /^^ ^Hydf. The probability to overcome a proposed MCMC 
jump between these two likelihood values determines the convergence of a MCMC sam- 
pler. The acceptance probability is therefore related to the difference of the log-likelihoods 
logLHt — logLHji. We want the hottest temperature to allow jumps within the posterior 
surface and we want this to happen about every, say, 1000 iterations. This number allows oc- 
casional jumps at the hottest temperature (smallest inverse temperature) which is therefore 
chosen to be 

_ log(lOOO) .^2) 
logLHt- log LH,- ^''^ 

We then use an exponential temperature scheme for k = 10 chains: 

/?.=/5M,J = {l,---,fc}- (13) 

For each iteration and each chain, new parameter values are proposed. Of course this is 
only meaningful when the current state of the chain is not the null-model. If the present 
state of the sampler is model Aig (or Ait depending on the comparison), independent normal 
distributions are chosen to propose new jumps. Pilot runs are first used to find appropriate 



proposal variances. The acceptance probability 



'or a proposed candidate is derived by 



35l | . The proposed swaps between arbitrary 



computing the Metropolis-Hastings ratio [15, 
pairs of chains are done in the way 
above. 

The transdimensional jumps between null-model Ain and model Aig (or Ait) are con- 
ducted by RJMCMC steps [271] ■ We implemented a death and a birth proposal which either 
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attempts to jump from model A4s (or Ait) to model or from model Ain to model A4s 
(or Mt). 



B. Reversible Jump Markov chain Monte Carlo 

The reversible jump approach requires a random variable r with distribution g(T) that 
matches the dimensions of the parameter space across models. In addition, a function is 
defined that does the dimension matching. In the present case it is a function based on 
death and birth events. The one-to-one transformation in the 'birth' transition creates a 
new signal with parameter vector 6' and has the form tQ^i{T) = t = 0' . The inverse 'death' 
transformation that annihilates the signal, tgii •= 'tit^o^ form ti^o(6') = 0' = t. The 
Jacobian of both transformations is equal to 1. The acceptance probability for the creation 
process is therefore 



where 6' is the parameter vector of the current existing signal. Since we chose equal prior 
probabilities for both models, we have p{M.s)/p{A4n) = p{A4n)/p{A4s) = 1- 

As one can see in Eq{TH for the creation process, the prior distribution p{0') at an existing 
parameter 6' is found in the enumerator while the proposal distribution q{0') at the existing 
parameter is present in the denominator. On the other hand, in EqJTSj for the annihilation 
process, the proposal value q{0) of a new proposed parameter vector 0' = t is found in the 
enumerator while the prior p{6') is in the denominator. This means that a larger parameter 
space, which naturally yields smaller prior values, results in more likely accepted deaths than 
accepted births. As a consequence, the sampler will prefer sampling from the null-model 
which means that the evidence of a signal will be smaller if we increase the parameter space. 
This makes perfect sense as we expect the evidence of a signal to fade if we integrate over 
a larger parameter space. 

The other fact that the proposal distribution enters on opposite sides of the fraction in 
EqJT^ and EqlTS] reveals the difficulty on the choice of the proposal distribution. In order 




(14) 



where 6' = t is drawn from ^(t). The annihilation process is in turn given by 




(15) 
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to understand the effect of the proposal distribution, let us consider the following three 
scenarios 

1. Suppose we do not know the major posterior mode and therefore choose the parameter 
vector of a new signal to be drawn from a wide spread proposal distribution. We could 
choose the proposal distribution to be the same as the prior distribution, in which case 
prior and proposal would cancel out. Prom the sampling point of view, the samples 
would account for the prior distribution and the acceptance probability would merely 
contain the likelihood. However, it would be unlikely to find the narrow mode in the 
likelihood surface by ineptly poking around in the entire parameter space restricted 
by the prior. Such a RJMCMC sampler would rarely accept jumps between the mod- 
els. Only very long runs would give sufficient information about what proportion the 
sampler naturally stays in which model in order to draw reliable conclusions. 

2. Suppose we wrongly assume the posterior mode to be concentrated in some area of 
the parameter space far away from the actual posterior mode. We would choose 
the proposal distribution to have the major probability in some wrong area of the 
parameter space. Naturally the draws from such distribution would privilege proposals 
in the wrong area of the parameter space but the acceptance probability would repress 
births and support deaths in the area as the proposal values would be naturally high. 
The sampler is balanced but mixing would be poor as proposes in the correct area of 
the parameter space are rare, even though their acceptance would be facilitated. The 
mixing would be even worse than in the first scenario and even longer runs would be 
needed to reveal reliable information about how long the sampler naturally stays on 
average in which model. 

3. Suppose we have a vague idea about where the major posterior mode is located and 
choose our proposal density to be centered around that area. The samples would be 
drawn preferential in that particular area of parameter space but the ratio of proposal 
and prior would compensate for that in the acceptance probability. However, the 
likelihood ratio in this area would have a major impact and the sampler is more likely 
to jump between the models revealing the proper ratio and model probability in a 
much shorter sampling period. 



13 



We see that the choice of a proper proposal distribution is very important. A good 
proposal distribution ought to have the major probability mass concentrated around the 
expected posterior mode but with long tails in order to cover the entire prior. This is 
achieved by a mixture distribution between a normal distribution with small variance and 
a uniform distribution that covers the prior range. If we are to compare the null-model Ain 
and the true model Ait, we are in a lucky position. The mean of the proposal distribution 
is most likely to be found around the true parameter values and we only need to find a 
variance in the same order of magnitude as the posterior mode which can be determined by 
pilot runs. 

Things are different when we compare null-model Ain and the simplified model Aig- The 
posterior mode can not be expected at the true parameter values as the wave form of the 
simplified model is definitely not best fit at the true parameter values. Just augmenting the 
variance of a proposal distribution with mean at the true parameter values would result in 
bad mixing. We therefore need pilot runs at higher signal-to-noise ratio in order to determine 
the vague center of the posterior mode in the simple model case. The information acquired 
from such runs serves to determine a suitable proposal distribution. 



C. Within-model Metropolis-Hastings sampling and re-parameterization 



iov chain 



35|. The 



The sampling process of the individual chains when the current state of the Mart 
is in the model that postulates a signal, is done by a common MH step 15|, 
proposals here is tailored to the expected posterior shape by choosing a very heavy tailed 
distribution. This was accomplished by mixing a normal distribution with exponentially 



varying variance 32|, |36|, l37|] . 



The high correlation of the mass parameters mi, and m2 in the posterior distribution 



needed to be accounted for 
1.5 PN time to coalescence 



jy re-expressing them in terms of the following Newtonian and 
38| 

8/3 {mi + m2f 



Ai = Fi ■ (mi + 1712) 



A2 = -F2 ■ (mi + 1712) 



-5/3 



mim2 
(mi + m2)' 

17111712 



(16) 



(17) 



where Fi = ^{t^/o) and F2 = f (tt/o) This turned out to work very well as 

re-parameterization technique for our sampler. 
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Since mi and 1712 can be written in terms of Ai and A2 according to the following expres- 
sions 



m, = ^{c^- sjd - ACl/' ) , (18) 



m, = \(c, + sjcl - ACl" ) , (19) 



where Ci = and C2 = |r j , it follows that the Jacobian of this transformation i 
equal to 



det J = '-^L^ ^ (20) 

(FiA2)2-4(F2Ai)2C2'/' 

Since in the original parameter space we defined a joint density for 
{mi, 7712, r}, the new joint prior distribution of Ai, A2, and r is given by 



p(Ai, A2,r)=< 



p(mi(Ai,A2),m2(Ai,A2),r)| det J\ if m„j„<mi(Ai, X2)<m^ax 

and m„i„<m2(Ai, A2) < (21) 
otherwise 



where mi(Ai, A2) and m2(Ai, A2) are given by Eq. [T51 and Eq. [13 

The sampling techniques described in the previous subsections are then used to sample 
from this new multidimensional parameter space {0, (Ai, A2, ^c, ^c, ''^Y^}- 

IV. DESCRIPTION OF THE SIMULATIONS 

For our simulations we have created data sets from ^^true" wave forms of three hypotheti- 
cal binary inspiral systems (BI). We will consider two scenarios where the "irue" wave form 
is either of PN 2.0 or PN 3.5 order. The detection of each scenario is attempted by either a 
PN 1.5 or a PN 3.0 wave form respectively. 

In the first scenario, PN3.0/3.5, "trae" model Ait containing a PN3.5 wave form is used 
for creating the observed data d{0; f) = St{6; f) + n{f). In the frequency domain a PN3.5 



signal has the form 
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39| 



~stie; f) = Af-'" exp [t {G{e- f) + H{e- /)^3.5(^; /))] (22) 

where 

f) = «i.5(0; /) + a2.o(^; /) + «2.5(0; /) + «3.o(0; /) + «3.5(e; /) (23) 
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with 



«i.5(^;/) 
«2.o(^;/) 

a2.5(^;/) 



1 + 



20 



9 V336 



_ / 3058673 



4M^ 
5429/x 



TT 



1016064 1008M 144M2 ) 
38645 38645 



756 
65/i 



log ( V6(7rM/) 



«3.o(^;/) 



9M 

11583231236531 
4694215680 



M 



1/3 

252 

3 log (v^(7rM/)^/^))^ (vrM/)^/^ 

6407r2 6848- 0.57721 \ 
~ 21 J 

1760-11831 



1533559782 7 2 25571^ 



3048192 
76055/i2 127825/1^ 



12 

6848 



TX 



1728M2 1296M3 21 
/ 77096675 378515/i 74045/i2 



V 254016 



1512M 756M2 



3 9240 
log(4) (vrM/)^/^ 



+ 



12320-1987 
~9 3080 

{nMff 



and 



and 



H{e-f) 



128 



-5/3 



G{e-J) = 2nftc- 



TT 



/4 



(24) 



(25) 



(26) 



with coalescence time tc and coalescence phase (pc, involved masses mi and m2, total mass 
M = nil + reduced mass = mim2/M, and chirp mass Mc = {rn^rn^l MY/^ . 



The amplitude A is related to the intensity of gravitational wave [38|] and in the stationary 
phase approximation, A oc mI^^ /r, where r is the distance between detector and source. 
The model is determined by the five parameters = {mi, m2, r, tc, 0c}"^- The proportional 
factor depends on the relative orientation between detector and source, and we will assume it 
to be constant and equal to 1 since we are more generally interested in the broad assessment 
of the evidence of a signal. 

The signal used for detection is the 3.0 PN approximation, which will serve as the sim- 
plified model Aig- The PN3.0 approximation formulated in the frequency domain is given 
by the following expression 



5,(0; /) = Af-'^' exp [z {0(0; f) + H{d; /)^3.o(^; /))] 



(27) 
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where 



^3.o{e;f) = ai.5(0;/) + a2.o(0;/) + «2.5(0;/) + «3.o(0;/) 



(28) 



In the exact same way the second scenario involving PNl. 5/2.0 wave forms is approached 
where the data are created using a supposedly "^rue" signal 



~st{e; f) = Af-'" exp [i {G{d; /) + H{e; /)V^2.o(0; /))] 



(29) 



where 



V'2.o(0;/) = ai.5(0;/) + «2.o(0;/). 



The simphfied model Ms uses the lower 1.5 PN expansion 



(30) 



~st{e; f) = Af-'/' exp [i (G(0; /) + H{e; f)a,,,{e; /))] 



(31) 



The distance of each binary system is varied in order to obtain different signal-to-noise 
ratios. The noise realizations are drawn in such a way that they correspond to the approx- 
imated expression of the one-sided power spectral density of initial LIGO citeChronopou- 
los:2001 

(32) 



Snif) = f 



7'- 



1+ I ^ 
JO 



with 5*0 = 8.0 X 10~^^ Hz~^ being the minimum noise of the detector and /o = 175 Hz the 
frequency at which the sensitivity of the detector reaches its maximum. 

The noise samples arc generated directly from the noise spectrum in the following way. 
Let us assume the noise to be white, Gaussian distributed N(0, 1) (i.i.d. standard normal at 
time t). Its finite [T^jTe] Fourier transform is given by 

n{f) = T{n{t)} = [ n{t)exp{-2mft)dt , (33) 



which we can write as 

Hn{t)} 



n{t) cos(27r/t)dt-i / n{t) sin(27r/t)dt 



:= R{f) - il{f) 



(34) 
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We aim for deriving expected values and variance of the real part R{f) = 
J^" n{t) cos(27r/t)dt and imaginary part /(/) = J^" n{t) sin(27r/t)dt of the Fourier trans- 
form. Since n{t) is normally distributed with zero mean, it follows that also the expectation 
values E = and E (/(/)) = 0. From this consideration it follows that the variances 

of the real and imaginary parts of the Fourier transform of the noise are equal to 



i.i.d. 




n{t) sm{2Tcft)dt 



n{t')n{t") sin(27r/t') sin(27r/t")dt'dt" 



\t) sm\2nft)dt] 



E{n^{t))sm\27rft)dt 



sm'^{2nft)dt 



and in the same way 



E(i?^(/)) 



cos\2nft)dt 



(35) 



(36) 



Ts 



while the expectation value of the product between the real and imaginary part of the noise 
is equal to zero. For large AT := - T, it is E[P{f)] ^ AT/2 and E[R^{f)] ^ AT/2. 
Therefore, the samples of n{f) can be generated by sampling n^-elf) ~ N(0,AT/2) and 
^im(/) ~ N(0,AT/2) according to a white spectrum of Gaussian noise of variance 1/2. 
These considerations indicate that we can generate the noise samples directly in the Fourier 
domain by sampling the real and imaginary parts of the noise from two independent random 
number generators that are Gaussian distributed and have both equal variance Sh{f)/2 
{Sh{f) being the one-sided power spectral density of the noise). 

Since the time required to perform a single-signal simulation was several days on the Jet 
Propulsion Laboratory Dell XEON cluster (running 1024 Intel Pentium 4 Xeon processors), 
we decided to perform only three simulations for three different binary systems. These were 
selected to be close to the "corners" of the (Ai, A2) region of the mass-space, corresponding 
respectively to the mass-pairs given in Tab. [B The reason why we chose binary systems 
of such particular mass constellations can be seen in the following Fig. [2] where the masses 
are drawn in the re-parameterized (Ai, A2)-plane. Here, the chirp mass is defined as = 
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system 


M, 


rj 


mi 


m2 


tc 


4>c 


r 


Bl 


O.879O5M0 


0.1875 


I.8M0 


O.6M0 


42.76933s 


0.2 rad 


16 


- 24Mpc 


B2 


3.O895O6M0 


0.0112931 


45. OM© 


O.52M0 


5.26426s 


0.2 rad 


22 


- 35Mpc 


B3 


31.85576M0 


0.24 


45. OM© 


30.0Mq 


0.10778s 


0.2 rad 


100 


- 250Mpc 



TABLE I: Table of the parameters of the three example binary systems B1-B3. 
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FIG. 2: Three example binary systems in the Ai, A2 plane. 

/M)i/5 and the mass function rj = ? which is the ratio between the reduced 

mass and the total mass of the binary system. Fig. [2] shows the three mass constellations close 
to the three corners of the (Ai, A2) triangular plane. The bottom left corner corresponds to 
large masses with low mass ratio, in the top left corner mass constellations are found with 
large and small masses (large mass ratio). Towards the right corner, the masses become 
small. 

V. RESULTS 

We created data sets for the three different binary systems given in Tab. [H and changed 
their distances in such a way that the resulting SNRs would give a detection probability 
varying within its extremes. This resulted in varying the SNR within the interval (3, 12). 
For the binary system Bl, we simulated 12 different distances of varying step width between 
16 and 24 Mpc. We found this step width to be sufficient in order to capture the variability 
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of the calculated detection probability as a function of the SNR. For binary system B2 
we similarly took 12 different distances between 22 and 35 Mpc, while for system B3 we 
considered 16 different distances in the range of 100 — 250 Mpc in steps of 10 Mpc. For 
each of the 40 distances considered we generated 20 different noise realizations, resulting 
in a total of 800 data sets. Our MCMC sampler was applied four times on each data set 
for covering the specified model comparisons. This yields a total of 3200 x 10 simulated 
Metropolis- Coupled MCMC chains, each of which was stopped after 300 000 iterations after 
inspecting that such a number was sufficient for our purpose. 

The simulated data were sampled at 4096 Hz for a duration of about 24 s, and they were 
produced by embedding the different signals into noise samples that where generated in the 
Fourier domain as described in Sec JIVI Since all MCMC runs were conducted after pilot 
runs at higher SNRs, the burn-in period was kept very short as the proposal distributions 
were optimized to the target distribution and mixing was very efficient. From the MCMC 
output we discarded just the first 10 000 iterations as burn-in, while short-term correlations 
in the chain were eliminated by "thinning" the remaining terms: every 100**^ item was kept 
in the chain. 

The integration bandwidth for the likelihood was chosen from 12 Hz up to the frequency 
of the last stable orbit or 600 Hz, whichever is the smaller. Since the SNR is negligible above 
600 Hz, we fixed this to be the upper frequency cut-off. For the Bl system, we derived a 
frequency of 1832.2 Hz at the last stable orbit which gives an integration limit of 600 Hz for 
Bl. This translates in 14355 complex samples that contribute to the posterior distribution. 
In the case of system B2 instead, the last stable orbit is at 96.6 Hz, resulting in 2065 complex 
samples involved in the determination of the likelihood function. Finally, for the high mass 
binary system B3 in our set of systems the last stable orbit is at 58.6 Hz, implying now only 
1138 complex samples over which the likelihood is calculated. 

After we conducted the MCMC simulations we derived the posterior detection proba- 
bilities for the competing models from the MCMC outputs regarding the three example 
binary systems, the four different model comparisons, and the different sets of SNRs. For 
each binary system we computed the posterior probabilities for the considered scenarios and 
contrasted the probability of detection based on a lower order PN expansion against the one 
using the true wave form. This is displayed in Figs. O HI [5] respectively. 

In order to derive the detection curves we computed the proportions of the states in which 
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the Markov chain recurred to the null model or the model containing a signal. This was 
done for the entire set of noise realizations for a given SNR. Lines connect the estimates of 
the posterior detection probabilities resulting in an interpolated function of the SNR. This 
is a monotonically increasing function of the SNR, reaching asymptotically 1 as the SNR 
goes to infinity. A common feature to these figures is the uncertainties due to the noise that 
the signal detection probability shows at a given SNR. In particular, these uncertainties are 
more pronounced when the gradient of the detection probabilities is at its maximum. Note 
also that the difference between the detection probabilities associated to the "true model" 
and the approximated one is much smaller than these uncertainties. The uncertainties are 
displayed as vertical bars: the 50% quartiles (thick bars), and the outer quartiles (thin lines) 
associated with the 20 noise realizations. The inner 50% quartiles are divided by a small 
line which represents the median. It is interesting to see in Fig. [3] that, in some cases, for a 
given binary system and SNR, detection probabilities as low as or as high as 1 are possible, 
merely on the effect of the noise realization. 



It is worth mentioning that analyses performed within the frequentist framework 
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38, 



401] and aimed at comparing the detectability of a signal by using a simplified wave form were 
focused entirely on estimating the resulting loss of SNR. The Bayesian model comparison 
presented here has the inherent ability to estimate probabilities and their uncertainties due 
to noise, providing much more insights into this issue. 

Another interesting feature shown by the detection probability curves is their asymptotic 
dependence on the SNR. While the probabihty of detection always converges to 1 as the SNR 
goes to infinity, it does not necessarily always goes to zero with the SNR. The reason for 
this lies in the Bayesian approach in which we assumed equal prior probability for both, the 
null-model and the model that contains a signal. In the Bayesian context, all probabilities 
represent a degree of belief. They are based on the prior information and on the information 
that is given by the data by means of the likelihood. The more data we have, the more 
new information we obtain from the posterior distribution. The less data we have, the more 
impact the prior has on the posterior. In an extreme scenario with no data at all, the 
posterior is equal to the prior. We used this fact to test the correctness of our RJMCMC 
sampler as it must sample properly from the prior distributions and prior model probabilities 
when no data are present. The three binary systems considered in this paper each imply 
different likelihoods. For example Bl, with its small masses, has a spectrum that nicely 
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falls into the part of the observable band of the detector where the instrumental noise is at 
its minimum. On the other hand, the system B3 shows an energy spectrum whose upper 
frequency cut-off is equal to 58.6 Hz, with a resulting 1138 frequency bins over which the 
likelihood is calculated. 

The diverse data sets are reflected in Figs. [31 HJ O For Bl with its 14355 involved samples, 
the detection probability converges to almost zero for low SNRs as the data provide sufficient 
evidence for the non-existence of a signal even though the prior suggests a probability of 
0.5. It is in the nature of the Bayesian approach that the scarcer data for B2 and B3 
provide less evidence resulting in a posterior detection probability of around 0.2 and 0.4, 
respectively, when the SNR approaches zero. A different choice for the prior probability 
on P{A4n) = 1 — P{-Ms) would change the course of the posterior detection curves but 
with increasing number of data samples and SNR, the likelihood dominates the posterior 
distribution. 

In order to point this up, we created a graph comparable to Fig. [5] (bottom) in which the 
results of B3 based on the PN3.0/3.5 comparison are shown with a different prior probability 
of P{Ain) = 0.99 for the null model. The model that contains a signal has consequently 
a prior probability of 0.01. The result is illustrated in Fig. [61 This plot corresponds to 
bottom graph of Fig. [5] with the only difference that a more pessimistic prior probability 
P{Ain) = 0.99 on the null model has been applied. When comparing Fig. [6] to Fig. [5] 
(bottom), we see that the detection probability is significantly lower at SNR < 7 due 
to the higher prior probability on the null model. However, at an SNR of around 6, the 
detection curve jumps up quickly in Fig.[6]and a detection probability of 1 is reached in both 
figures roughly at an SNR of 8 because the evidence of a signal in the data is overruling the 
prior probability in both cases. 

We will now focus on the parameter estimates. We have compiled plots which address 
the impact of the use of lower PN order wave forms on the bias of posterior distributions 
of the parameters. Along the lines of Figs. [31 [H, [5], we plot the posterior distributions of 
the parameters when the posterior detection probability reaches a value of 1. The posterior 
distribution of the parameters is hereby an integration over the noise by incorporating all 
20 realizations of the MCMC outputs. Each output corresponds to the SNRs at which the 
detection probabilities in Figs. [31 [2 [SI reach their maximum. We only concentrated on the 
four parameters mi, m2, r, tc as the phase (pc is of no particular physical interest. We 
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displayed the posterior density of the masses as a 2D joint probabihty density in the form 
of a contour plot. We computed the two-dimensional 50% and 95% credibility regions. 

For the three considered binary systems we generated a total of 36 plots for the distribu- 
tions of the chirp mass and the mass function i], as well as density plots for the distance 
r and the time to coalescence tc- We have chosen to plot the joint posterior probability 
of the mass parameters in the (Mc, ?7)-space because they are not as much correlated as 
(mi, 7712) in their posterior, which produce hard to visualize posterior densities. Although 
we could plot the posterior in the (Ai, A2) space, the joint posterior probability of Mc and r] 
is physically more meaningful. We divided the 36 plots into three sets corresponding to Bl 
(Fig. ED, B2 (Fig. ED, and B3 (Fig. ED . 

Figs. El [HI El display the true parameter values chosen in our simulations (dashed line in 
the case of r and tc-, and intersection of two dashed lines in the Mc,?7 plots). From visual 
inspection we notice that when the model matches the signal present in the data the posteri- 
ors cover well the true parameter values. In Fig. El however, the joint posterior distribution 
of the mass parameters are offset from the true values in the PNl. 5/2.0 comparison. The 
PN3.0/3.5 detection, on the other hand, reveals a much smaller offset for this particular 
signal. However, this is not true in general, as it can be seen for the PN3.0/3.5 comparison 
shown in Fig. [HI The offsets of the posterior distributions from the true values of the mass 
parameters are very obvious in both, the PNl. 5/2.0 and PN3.0/3.5 comparisons. The pos- 
terior is shifted over several of its standard deviation. Very striking is also the error in the 
time to coalescence for the PN3.0/3.5 comparison in the B2 signal. The mass parameters 
Mc and t] and time to coalescence are obviously the parameters subject to biases when using 
a simplified model. This is physically understandable since these three parameters define 
the phase of the signal. The posterior distributions of the mass-related parameters shown in 
Fig. [9] reveal smaller offsets. This is because the spread of the posterior distribution is less 
pronounced and the bias is therefore smaller compared to the posterior standard deviation. 
The posterior distribution of time to coalescence, however, is strongly offset from the true 
parameter value in the PN3.0/3.5 comparison. 

The graphical output only serves as a visualization. For an honest comparison, numbers 
are needed. To this end, in Tab. [Til "we show the true values, the 95% posterior credibility 
interval, the median and the mean of the chirp mass M^, based on the MCMC outputs. 
Tab. mil Tab. IIVI and Tab. [V] show the same entries for the parameters 77, tc-, and r. 
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simulation 
identification 


true value 


95% credibility 
interval (CI) 


posterior 
mean 


posterior 
median 


true value falls 
into 95% CI 


Bl: PNl.5/2.0 


0.87905 


[0.87867, 0.87978] 


0.87908 


0.87911 


/ 


Bl: PN2.0/2.0 


0.87905 


[0.87865, 0.87955] 


0.87906 


0.87907 


/ 


Bl: PN3.0/3.5 


0.87905 


[0.87855, 0.87942] 


0.87900 


0.87899 


/ 


Bl: PN3.5/3.5 


0.87905 


[0.87863,0.87942] 


0.87905 


0.87904 


/ 


B2: PNl.5/2.0 

/ 


3.08951 


[3.08066,3.11047] 


3.09384 


3.09576 


/ 


B2: PN2.0/2.0 


3.08951 


[3.07414,3.10998] 


3.08931 


3.09047 


/ 


B2: PN3.0/3.5 


3.08951 


[3.06665, 3.09953] 


3.07790 


3.08828 


/ 


B2: PN3.5/3.5 


3.08951 


[3.07999, 3.09852] 


3.08955 


3.08969 


/ 


B3: PNl.5/2.0 


31.85576 


[29.54749, 32.09530] 


31.01141 


30.97721 


/ 


B3: PN2.0/2.0 


31.85576 


[30.76357, 33.50058] 


31.90688 


31.95160 


/ 


B3: PN3.0/3.5 


31.85576 


[31.09861,34.68858] 


32.53342 


32.60697 


/ 


B3: PN3.5/3.5 


31.85576 


[30.52275,33.82286] 


32.08128 


32.07732 


/ 



TABLE II: Simulation results of chirp mass Mc- 



respectively. The right-most column of the tables compares whether the true values of the 
binary systems fall into the corresponding 95% posterior credibility intervals. 

The results seen in these tables show that the mass function, 77, and the time to co- 
alescence, tc, are the parameters that are most biased when estimated with a simpli- 
fied model. Note, first of all, that in all cases, the true values fall into the 95% cred- 
ibility intervals when the estimation is based on the true model. For the distance r 
and the chirp mass Mc the 95% credibility intervals cover the true values in all compar- 
isons. However, when applying a simplified model, the 95% credibility interval of the mass 
function 77 does not cover the true value in three cases: {BlrPNl. 5/2.0, B2:PN1. 5/2.0, 
and B2:PN3.0/3.5)}. In the case of the time to coalescence tc, we find instead four 
cases: {B2:PN1.5/2.0,B2:PN3.0/3.5,B3:PN1.5/2.0,B3:PN3.0/3.5}. Combining these cases, 
we have a total of 5 out of the 6 simple model comparisons (PNl.5/2.0 and PN3.0/3.5) that 
fail to retrieve all their parameters within the 95% credibility region. The only simple model 
comparison that yields 95% credibility intervals that overlap all the true parameter values 
is Bl: PNl.5/2.0, although this is only marginal. 
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simulation 
identification 


true value 


95% credibility 
interval (CI) 


posterior 
mean 


posterior 
median 


true value falls 
into 95% CI 


Bl: PNl.5/2.0 


0.18750 


fO.21344, 0.244121 


0.22521 


0.22591 




Bl: PN2. 0/2.0 


0.18750 


[0.17844,0.20027] 


0.18783 


0.18819 


/ 


Bl: PN3.0/3.5 


0.18750 


[0.18110,0.19046] 


0.18559 


0.18564 


/ 


Bl: PN3.5/3.5 


0.18750 


[0.18292,0.19184] 


0.18758 


0.18753 


/ 


B2: PNl.5/2.0 


0.01129 


[0.01371,0.01525] 


0.01435 


0.01449 




B2: PN2.0/2.0 


0.01129 


[0.01070,0.01210] 


0.01129 


0.01133 


/ 


B2: PN3.0/3.5 


0.01129 


[0.01077,0.01101] 


0.01086 


0.01095 




B2: PN3.5/3.5 


0.01129 


[0.01121,0.01137] 


0.01129 


0.01129 


/ 


B3: PNl.5/2.0 


0.24000 


[0.22428, 0.249981 

L ' J 


0.24350 


0.24118 


/ 


B3: PN2.0/2.0 


0.24000 


[0.22821,0.24998] 


0.24436 


0.24245 


/ 


B3: PN3.0/3.5 


0.24000 


[0.22796, 0.24999] 


0.24445 


0.24250 


/ 


B3: PN3.5/3.5 


0.24000 


[0.22590, 0.24999] 


0.24426 


0.24200 


/ 



TABLE III: Simulation results of mass ratio rj. 



simulation 
identification 


true value 


95% credibility 
interval (CI) 


posterior 
mean 


posterior 
median 


true value falls 
into 95% CI 


Bl: PNl.5/2.0 


42.76933 


[42.76716,42.77022] 


42.76854 


42.76856 


/ 


Bl: PN2.0/2.0 


42.76933 


[42.76805,42.77089] 


42.76937 


42.76940 


/ 


Bl: PN3.0/3.5 


42.76933 


[42.76931,42.77313] 


42.77113 


42.77115 


/ 


Bl: PN3.5/3.5 


42.76933 


[42.76740,42.77111] 


42.76937 


42.76935 


/ 


B2: PNl.5/2.0 


5.26426 


[5.23179, 5.26308] 


5.24587 


5.24686 




B2: PN2.0/2.0 


5.26426 


[5.24868, 5.28231] 


5.26398 


5.26439 


/ 


B2: PN3.0/3.5 


5.26426 


[5.52561,5.58354] 


5.54898 


5.54958 




B2: PN3.5/3.5 


5.26426 


[5.24267,5.28606] 


5.26438 


5.26467 


/ 


B3: PNl.5/2.0 


0.10778 


[0.08805,0.10761] 


0.09706 


0.09723 




B3: PN2.0/2.0 


0.10778 


[0.09796,0.11714] 


0.10802 


0.10794 


/ 


B3: PN3.0/3.5 


0.10778 


[0.13103,0.15565] 


0.14444 


0.14390 




B3: PN3.5/3.5 


0.10778 


[0.09220,0.12211] 


0.10999 


0.10884 


/ 



TABLE IV: Simulation results of time to coalescence tc- 
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simulation 
identification 


true value 


95% credibility 
interval (CI) 


posterior 
mean 


posterior 
median 


true value falls 
into 95% CI 


Bl: PNl. 5/2.0 


16.00 


[14.04,21.631 


17.00 


17.22 


/ 


Bl: PN2.0/2.0 


16.00 


[14.03,21.55] 


16.96 


17.17 


/ 


Bl: PN3.0/3.5 


16.00 


[14.18,22.04] 


17.20 


17.44 


/ 


Bl: PN3.5/3.5 


16.00 


[14.15,21.92] 


17.18 


17.40 


/ 


B2: PNl. 5/2.0 


22.00 


[19.13,31.17] 

L ' J 


23.67 


24.07 


/ 


B2: PN2.0/2.0 


22.00 


[19.02,30.98] 


23.61 


23.98 


/ 


B2: PN3.0/3.5 


22.00 


[19.23,31.83] 


23.86 


24.59 


/ 


B2: PN3.5/3.5 


22.00 


[19.23,31.35] 


23.71 


24.11 


/ 


B3: PNl. 5/2.0 


100.00 


[84.54, 168.42] 


111.55 


115.45 


/ 


B3: PN2.0/2.0 


100.00 


[86.07, 174.08] 


114.05 


118.16 


/ 


B3: PN3.0/3.5 


100.00 


[86.94, 172.77] 


114.94 


118.91 


/ 


B3: PN3.5/3.5 


100.00 


[85.92, 169.58] 


113.05 


116.76 


/ 



TABLE V: Simulation results of distance r. 



In summary, we see that the bias in the estimated parameters based on a simpler model 
is larger than the statistical uncertainty. However, we should note that the SNR we have 
been considering corresponds to the value at which the posterior detection probability just 
reaches the value of one. Since the statistical error is a monotonically decreasing function 
of the SNR while the bias is not, we conclude that the difference between statistical and 
systematic error increases for larger SNRs. 

These results reveal that parameter estimates based on simplified models are not very 
reliable, since the systematic error is higher than the uncertainty of the posterior distribution. 
Furthermore, the use of higher order post Newtonian wave forms does not abate this problem, 
as it has been shown in Fig. [HI 

VI. CONCLUSION 

We have shown that, within the Bayesian framework, the probability of detection is not 
impeded by using a simplified model for detecting wave forms of higher PN order in the low- 
SNR regime. The Bayesian approach provides the means to gain insight into the variation of 
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the detection probability over different noise realizations. We have shown that the difference 
between the posterior detection probabilities corresponding to the true and the simplified 
model is very small as compared to its variance over different noise realizations. We have 
further shown that the systematic error in the Bayesian estimates, on the other hand, can 
be larger than the statistical uncertainties. This is also in agreement with results obtained 
within the frequentist approach, discussed in the literature by others 
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38|,|40(. However, 



it is based on finding the best fit of the involved wave forms while in our Bayesian framework 
an integration is performed over the entire posterior distribution which implies detection and 
estimation simultaneously. We can therefore analyze the posterior distributions, conditioned 
on the model that involves a signal which provides us with credible estimates in the low-SNR 
regime. 

We find that the estimates of rj and tc based on simplified models need to be taken with 
caution in both the PNl.5/2.0 and in the PN3.0/3.5 the offset is unpredictable. 

The only parameter that is accurately recovered throughout our simulations is the distance 
r which is clear as it only appears in the amplitude term and is not affecting the phase 
evolution of the signal. The chirp mass could also be retrieved within the 95% credibility 
intervals but yet shows a visible offset. With increasing SNR, however, the statistical error 
becomes smaller while the systematical offset remains constant. Given these findings we 
conclude that post Newtonian approximations, regardless of order, can be precarious for 
detecting "true" gravitational wave forms. 
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FIG. 3: Probability detection curves of Bl for the PNl. 5/2.0 and PN2. 0/2.0 comparisons (a) and 
the PN3.0/3.5 and PN3.5/3.5 comparisons (b). The vertical gray bars indicate the 50% quartiles, 
and the thin lines refer to the outer quartiles associated with the 20 noise realizations. The inner 
50% quartiles are divided by a small line which represents the median. The lower PN vs. higher 
PN comparisons are shown as solid lines (detection curves) and light gray bars (quartiles) while 
the equal PN comparisons are displayed as da^id lines and dark gray quartile bars. 
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FIG. 4: Probability detection curves of B2 for the PNl. 5/2.0 and PN2. 0/2.0 comparisons (a) and 
the PN3.0/3.5 and PN3.5/3.5 comparisons (b). The vertical gray bars indicate the 50% quartiles, 
and the thin lines refer to the outer quartiles associated with the 20 noise realizations. The inner 
50% quartiles are divided by a small line which represents the median. The lower PN vs. higher 
PN comparisons are shown as solid lines (detection curves) and light gray bars (quartiles) while 
the equal PN comparisons are displayed as da^^d lines and dark gray quartile bars. 





FIG. 5: Probability detection curves of B3 for the PNl. 5/2.0 and PN2. 0/2.0 comparisons (a) and 
the PN3.0/3.5 and PN3.5/3.5 comparisons (b). The vertical gray bars indicate the 50% quartiles, 
and the thin lines refer to the outer quartiles associated with the 20 noise realizations. The inner 
50% quartiles are divided by a small line which represents the median. The lower PN vs. higher 
PN comparisons are shown as solid lines (detection curves) and light gray bars (quartiles) while 
the equal PN comparisons are displayed as da^^d lines and dark gray quartile bars. 
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FIG. 6: Probability detection curves of B3 for the PN3.0/3.5 and PN3.5/3.5 comparisons. The 
vertical gray bars indicate the 50% quartiles, and the thin lines refer to the outer quartiles associated 
with the 20 noise realizations. The inner 50% quartiles are divided by a small line which represents 
the median. The lower PN vs. higher PN comparisons are shown as solid lines (detection curves) 
and light gray bars (quartiles) while the equal PN comparisons are displayed as dashed lines and 
dark gray quartile bars. 
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FIG. 7: MCMC generated posterior densities for Bl. Part (a) shows the comparison with data of 
a PN2.0 wave form and (b) uses data with the PN3.5 wave form. Each of the figures in (a) and 
(b) show two different model comparisons based either on a lower PN order signal or the same PN 
order that was used in the data. The left column shows the joint posterior density of the mass 
parameters Mc and r/ in form of the 95% credibility area that contains 95% of the probability 
mass and the inner 50% credibility region colored in gray. The middle column shows the MCMC 
generated kernel density estimate (KDE) of the distance r and the right column the KDE for the 
time to coalescence tc- The true parameters values are indicated as dashed lines. 
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FIG. 8: MCMC generated posterior densities for B2. Part (a) shows the comparison with data of 
a PN2.0 wave form and (b) uses data with the PN3.5 wave form. Each of the figures in (a) and 
(b) show two different model comparisons based either on a lower PN order signal or the same PN 
order that was used in the data. The left column shows the joint posterior density of the mass 
parameters and r/ in form of the 95% credibility area that contains 95% of the probability 
mass and the inner 50% credibility region colored in gray. The middle column shows the MCMC 
generated kernel density estimate (KDE) of the distance r and the right column the KDE for the 
time to coalescence ic- The true parameters values are indicated as dashed lines. 
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FIG. 9: MCMC generated posterior densities for B3. Part (a) shows the comparison with data of 
a PN2.0 wave form and (b) uses data with the PN3.5 wave form. Each of the figures in (a) and 
(b) show two different model comparisons based either on a lower PN order signal or the same PN 
order that was used in the data. The left column shows the joint posterior density of the mass 
parameters Mc and r/ in form of the 95% credibility area that contains 95% of the probability 
mass and the inner 50% credibility region colored in gray. The middle column shows the MCMC 
generated kernel density estimate (KDE) of the distance r and the right column the KDE for the 
time to coalescence tc- The true parameters values are indicated as dashed lines. 
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